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Abstract 

This work originates from a heart's images tracking which is to 
generate an apparent continuous motion, observable through intensity 
variation from one starting image to an ending one both supposed 
segmented. Given two images po and pi, we calculate an evolution 
process p(t, •) which transports po to p\ by using the optical flow. In 
this paper we propose an algorithm based on a fixed point formulation 
and a space-time least squares formulation of the transport equation 
for computing a transport problem. Existence results are given for 
a transport problem with a minimum divergence for a dual norm or 
a weighted -ffg-semi norm, for the velocity. The proposed transport 
is compare with the transport introduced by Dacorogna-Moscr. The 
strategy is implemented in a 2D case and numerical results are pre- 
sented with a first order Lagrange finite element, showing the efficiency 
of the proposed strategy. 
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1 Introduction 

Modern medical imaging modalities can provide a great amount of infor- 
mation to study the human anatomy and physiological functions in both 
space and time. In cardiac magnetic resonance imaging (MRI) for example, 
several slices can be acquired to cover the heart in 3D and at a collection 
of discrete time samples over the cardiac cycle. From these partial obser- 
vations, the challenge is to extract the heart's dynamics from these input 
spatio-temporal data throughout the cardiac cycle [14], |16j . 
Image registration consists in estimating a transformation which insures 
the warping of one reference image onto another target image (supposed 
to present some similarity). Continuous transformations are privileged, the 
sequence of transformations during the estimation process is usually not 
much considered. Most important is the final resulting transformation and 
not the way one image will be transformed to the other. Here, we consider 
a reasonable interpolation process to continuously map the image intensity 
functions between two images in the context of cardiac motion estimation 
and modeling. 

The aim of this paper is to present, in the context of optical flow, an al- 
gorithm to compute a time dependent transportation plan without using 
lagrangian techniques. 

The paper is organized as follows. The introduction is ended, by recalling 
the optical flow model (OF) . In section [21 the algorithm is presented, and its 
convergence is discussed. In section [3] it is shown that the solutions obtained 
with the proposed algorithm are the solutions minimizing the same energy 
than the time dependent optimal mass transportation problem. Section U] 
is devoted to numerical results. In particular a 2D cardiac medical image is 
considered. 

1.1 The optical flow (OF) method 

Let p be the intensity function, and v be the velocity of the apparent motion 
of brightness pattern. An image sequence is considered via the gray-value 
map p : Q = (0,1) x Q — > M. where Q C M. d is a bounded regular domain, 
the support of images, for d = 1,2,3. If image points move according to the 
velocity field v : Q — > M. d , then gray values p(t,X(t,x)) are constant along 
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motion trajectories X(t,x). One obtains the optical flow equation. 



—p(t, X(t, x)) = d tP (t, X(t, x)) + (v\ V x p(t, X(t, x)) ) Rd = 0. (1) 



The previous equations lead to an ill-posed problem for the unknown (p, v). 
Variational formulations or relaxed minimizing problems for computing jointly 
(p, v) have been first proposed in [4] and after by many other authors. Here 
our concern is somewhat different. Finding (p,v) simultaneously is possible 
by solving a mass transport problem. Similarly to the work developed in 
[5j [6], a characterization of (p,v) as solution of a minimizing problem is 
developed. 

Let po an d pi be the cardiac images between two times arbitrary fixed to 
zero and one, the mathematical problem reads: find p the gray level function 
defined from Q with values in [0, 1] verifying 



The velocity function v, is determined in order to minimize the functional. 



Thus we get an image sequence through the gray-value map p. Let us 
mention [3], for example, where the optimal mass transportation approach 
is used in images processing. In this work, the optimal transport problem 
in 2D is decomposed in several ID optimal transport problems which are 
easier to numerically solved. The algorithm proposed here, is based on a 
least squares formulation for the transport equation ([9] for example), which 
differs from the methods proposed in [3] , or in [5] . Notice that the proposed 
transport in this paper differs from the optimal transportation, studied e.g. 
in the book of C. Villani [20]. 

2 Algorithm for solving the optical flow 

In this section, an algorithm is presented to solve the optical flow given by 
equations ([2]), and ([3]). Let us first specify our hypotheses. 

HI The domain SI is a bounded C 2,a domain satisfying the exterior sphere 
condition. 



d t p(t, x) + ( {v(t, x) | Vp(t, x) ) = 0, (i, x) e (0, 1) x 
p(0,x) = p (x); p(l,x) = pi{x)x e O 



(2) 




(3) 



3 



H2 The functions pi £ C 1,a (Q) for i = 0, 1, with po = Pi ° n 9CI. Moreover 
there exist two constants such that < j3 < pi < /3 in fi. 

Let p° £ C 1,a ([0, 1] x O) be given by p°(t,x) = (1 - t)po(^) + *pi(ac). We 
have ||^P ||co,«([o,i] x TT) ^ C(p ,pi) and 9tp°|an = 0. 

For each t £ [0,1], our need to solve problem ([2])-(j3]) is a velocity field 
vanishing on d£l. To do so, the following method is used. 
Assume that p n £ C 1,a ([0, 1] x $1) is given with p n (0, x) = po(x), and 
p n (l,x) = px(x) for all x £ 0. 

• Compute r/ n+1 such that 

-div(p n (V)V?f +1 ) = in n 

dn n+1 (4) 
P n (i,-)^ r = l ondfi, 



and set C n+1 (t) = — L / d t p n r) n+1 dx. 
\9n\ J n 

For each t £ [0, 1] compute <p n+1 solution of 



- div(p n (i, -)V(p n+1 ) = d t p n (t, •), in n 
ip n+1 = C n+1 (t) ondtt. 



(5) 



• Set v n+1 = Vip n+1 . 

• Compute p n+1 , L 2 -least squares solution of 

r 9 t p n+1 (t , x) + ( u n+1 (t,x)\ v P n+1 (t,x)) =0, (t, x) £ (o, i) x n, 

\ p n+l (0,x) = p (x); p n+l {l,x) = Pl {x) x £ ft. 

(6) 

Remark 2.1 If the requirement Po\qq = PiI^q ^ s canceled in hypothesis H2, 

dif n+l 

then the boundary condition of problem |5l) is replaced by — = and 

on 

we do not need anymore the constants C n+1 . 

For each t £ [0,1], since p n (t,-), and d t p n (t,-) £ C°' a (Tl), theorem 6.14 p. 
107 of [T3] applies, and there exists a unique (p n+1 (t, •) £ C 2,a (VL) solution of 
problem ([5]). In problem ([5]) the time t is a parameter. Since p n £ C 1,a , and 
dtf) n e ^0,0^ and C n+i e C ro,a j the c i assica i <7 2 ' a (ft) a priori estimates for 
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solutions to elliptic problems allow us to show that ip n+1 is a C°' a function 
with respect to time. So we have 

\\<p n+ llc°.«([o,l];C 2 >«(n)) < M (\\C n \\co^([o,i}) + II^P n |lc*°. a ([o,i]xIT))- 

Consider the extension of ip n+l by C n+1 outside of the domain f2, still de- 
noted byy? n+1 . Since the right hand side of equation ([5]) vanishes on dQ, 
this extension is regular, so the function v n+l vanish outside Vt and belongs 
to C°' a ([0,l];C 1 ' Q (R 2 )). 

Define the flow Xl +1 (s,t,x) £ C 1 ' a ([0, 1] x [0, 1] x M 2 ;]R 2 ) by 

f ^ +1 (s,t,x) = +^ +1 ( S ,^ +1 ( S ,t,x)) in (0,1) (?) 
1 x\ + \t,t,x) = x. 

This flows will be constant for x € d£l. Observe that solving the transport 
equation ([6]) in the least-squares sense is equivalent to solve this equation 
on each integral curves defined by (J7J. 

Set r(s) = p n+1 (s, X™ +1 (s, t, x)), and express equation (jfij) along the integral 
curves of equation (|7j). The equation is reduced to the following ordinary 
differential equation with initial and final conditions. 

(s*> = ° (8) 

I r(0) = Po(X!f +1 (0, t,x)); r(l) = f, i)). 

The L 2 least squares solution of (|SJ) minimizes -7-r(s), and is given by: 

r(a) = (1 - S )p (^ +1 (0, t, x)) + s Pl {Xl+\l, t, x)). 
Therefore the following representation formula for the function p n+1 is proved. 



Lemma 2.2 The 1? -least squares solution of problem ([6]) is given by 

p n+1 (t,x) = (l-t)p (Xl +1 (0,t,x))+t Pl (Xl +1 (l,t,x)). (9) 

Remark that the regularity of the function p n+1 is a consequence of the 
regularity of the flow X™ +1 . 

Let us now consider the convergence of the algorithm (UJ)-©. 
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Theorem 2.3 There exist (p,tp) G C 1 ([0,1] x x C°([0, 1]; C 2 (0)), 

L 2 -least squares solution, respectively solution of 



d t p+ ( Vy>| Vp) = 0, in(0,l)xO 
p(0,x) = po(a?)l = Pi(x) inO 

-div(p(t, -)V^) =5tp(t,-), infi 
ip = C(t) and V<p = on9fl 

mi/i C(i) defined as follow. 

-div(p(t, -)V??) = infi 
p(t, •) <9 n ?? = 1 on SO 

C = TanT / d tPVdx. 
\dSl\ Jn 



(10) 
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;i2) 



n n+l| 



Proo/. Since ||^ ||c .«([o,i])+II^P°llc7 .«([o,i]xn) is bounded, \\ip r ' ■ |ico,a([o,i];C 2 .° 
and |b n+1 ||c°' Q ([0,i]:C 1 ' a (R 2 )) are uniformly bounded in n. 
From lemma I2T21 there exists a unique p n+1 , the L 2 -least squares solution of 
(J6j) . Let us give an estimate for D 3 X™ +1 . Starting from 

D^l+^s, t, x)) = v n+1 (s, Xl +l (s, t,x)), 

we deduce (see [T]) 

D^X^is^x) = D 2 v n + 1 (s,Xl +1 (s,t,x))D 3 Xl +1 (s,t,x) 
D 3 X r ; +1 {t,t,x) = Id. { J 

Since D 3 D x Xl +1 {s,t,x) = P>iP> 3 ^+ +1 (s, *, ac) we get 

D 3 ^ +1 ( S ,t,rE) = e -//^(^+ 1 (r 1 X ; + 1 (r,t 1 x)))<ir /d _ (u) 

Thus ||D3U™ +1 ||cfo,a([o i] 2 xR 2 ) i s uniformly bounded in n. Moreover we have 

m 

D 2 Xl +l (s,t,x) = { v n+1 (s, t, x) | D 3 Xl +1 (s,t,x)) 

so ||-D2W n+1 ||c°' ce ([o,i] 2 x]R2) is bounded independently of n. 

From theorem 12.21 we deduce that ||p n ||ci,«([o i] x ?T) * s uniformly bounded. 

Since the embeddings 

C°' a ([0, 1]; C 2 ' a (Q)) C°([0, 1]; C 2 (ft)) and C^QO, 1] xTT) C^QO, 1] xO) 

are relatively compact there is a subsequence of (p n , ip n ) solution of (J3J)- dEJ) , 
still denoted by (p n , ip n ) converging to (p, tp) in C\[0, 1] xO) xC7°([0, 1]; C 2 (Tl)) 
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and (p, ip) is the solution of (|10p - (|12p provided the boundary conditions to 
be justified. The condition V(p n \gn = is valid for the approximations 
ip n (since the functions can be extended by C n outside of 0). So the con- 
vergence in C°([0, 1]; C 2 (Q)) yields the condition for the gradient of limit 
function. For the approximations of function p, the formula given in lemma 
12.21 combined with the regularity result show that the boundary conditions 
are exactly satisfied. These conditions are thus valid for the limit function 
due to the convergence in C 1 . □ We will show in the next section that 
the solution minimizes the same energy as for the time dependent optimal 
transportation mass. 



3 Interpretation of solutions to problem (flAJl) - f[T2T) 



In this section it is shown that the solution to problem (|10p -(|12p is a solution 
to a time dependent mass transportation problem. 

Zero is a bound from below of the following functional j Jq ||3tU+div(uV(^>— 
C))||^2(q) dt to be minimized with respect to (u,ip). 

1 f 1 

= i J ^ &tp + div(/9V((/9 " c ^ 2 L 2 m dt 

thus (p,ip — C) solution of (fTUj) - (fT2"j) minimizes 

1 I' 1 

Min - / \\d t u + dw(uX7i[j)\\ 2 L 2, n) dt. 

d t u+(u I W)=0 
«(0)=po; u(l)=pi in SI} 

So the solution (p,ip) of problem (fl0]l - (fT2l) . satisfies 

1 f 1 

(p,ip-C)= Argmin -/ \\d t u + div(uVip)\\ H -i dt. (15) 

{V>e£ 2 ((o,i);iri(fi)), Jo 

u€L 2 ((0,l);L 2 (U)), 
dtu+(u\V4>)=0, 
n(0)=po; u(l)=p\ in Q} 
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Theorem 3.1 Let (p,ip) be the solution of problem (|10p -(jl2 p given in the- 
orem [273\, then it satisfies 

(pj'Vip) = Argmin / / u\\v || 2 dxdt. (16) 

{«6L 3 ((o,i)i(H 1 (n)) a ) l Jo Jn 

MGL 2 ((0,l);i 2 (n)), 
dtu+div(uv)=0, 
d t u+(v\ Vu)=0, 
n(0)=po; u(l)=pi in Q} 

Proof. Observe that in problem (|16|) . the transport equation is solved with 
a L 2 -least square procedure. This is equivalent to find £ = u — (l — t)po + tpi 
such that 

dtC + (v | V£) = P R (pi - po + (v | V ((1 - t)po + t Pl ) ) 

where Pr is the L 2 projection onto the range of the transport operator for 
the velocity v. 

For u € L°°(fi), the expression v >— > div(uv) is well defined as a linear 
continuous on Hq(Q,). If moreover < (3 < u < (3 in Q, let H = Hq(Q) be 
equipped with the following inner product. 



(9,i/j)= [ u(V6\V^) dx, 
Jn 



which induces a semi-norm which is equivalent to the i/ 1 -norm. The Riesz 
theorem claims that for the linear continuous form 

£«C0) =< -diV {uv),Tp > H ;H>, 

there is a unique 9 G H such that 

C v {ip)= [ u ( ve | vv> ) dx, w € H. 
Jn 

Therefore v = V# for a 9 G -ffo(^) anc ^ problem ([16]) is reduced to 

(p,V(p)= Argmin / f u\\ V^|| 2 dxdt. (17) 

W>eL 2 ((0,i);Ho(^)rW 2 (n)), ^° ^ 
«6L 2 ((0,l);i 2 (n)), 
9tti+div(nV'0)=O, 
9tu+(V^|u)=0, 
n(0)=po; u(l)=pi in C} 
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or 



(p,Vip)= Argmin f [ u\\ VV'H 2 dxdt. (18) 



Since 



{-V>6L 2 ((0,l);^(Q)ntf 2 (n)), 

«eL 2 ((o,i);i 2 (n)), 

dtu— div(nV'0)=O, 
9 t u-(w| Vi/))=0, 
«(0)=po; w(l)=pi in 57} 

/ u||V^|| 2 cte = || -div(uV^)||^ 



2 

l ■ 



problem (fTTJ) reads: 



(p, V(f) = Argmin / || div(— uVi())\\ H -x dt. 



(19) 



{-^eL 2 ((0,i);H ( 5(Q)n// 2 (Q)), 
u£L 2 ((0,l);L 2 (Q)), 
dtu— div(uVi/>)=0, 
atu-(u|V^)=0, 
«(0)=po; u(l)=pi in SI} 



From the definition of the linear form £, observe that 

-|| div(-uV^) + d t u\\ 2 H -i = || div(-uV^)||#_i, 
so problem (|19p reads: 



1 Z" 1 

(p, V<p) = Argmin - / || div(iiV^) + d t u\\ 2 H -i dt. (20) 

{4>eL 2 {(o,i);Hi(n)nH 2 {n)), Jo 

nGL 2 ((0,l);i 2 (Sl)), 
8tu+div(uVtp)=0, 
d t u+{u\Vip)=0, 
u(0)=po; u(l)=pi in Q} 

Gathering (|15p with the previous result proves the theorem. □ 

Remark 3.2 The proposed transport is a one which minimizes the diver- 
gence of the velocity in a weighted dual norm of Hq 

Remark 3.3 The Dacorogna-Moser transport fT®/ p(t,x) = (1 — t)po(x) + 
tpi(x) with 

-Aip = d t p in ft; 
d n <p = 0; on dQ; 



9 



and the velocity v = ^ satisfi 



e.s 



(p,V(p) = Argmin f f \\v\\ 2 dxdt. (22) 

neL 2 ((o,i);i 2 (n)), 

9tu+div(u)=0, 
n(0)=po; w(l)=pi in f2} 

Indeed, as before, consider the space H = Hq equipped with the previous 
semi-norm, and set 

C v (tp) =< - div(v), ip > H . H ,=< d t U,1p >H;H' ■ 

Then the problem is reduced to 

1 f 1 

(p,Vy>)= Argmin - / \\d t u - div(V^)||^_i dt. (23) 

{V>ei 2 ((o,i);i? 1 (n)nJ3' 2 (n)), Jo 

u£L 2 {(0,l);L 2 {n)) 
a t u+div(VV')=0 

Using the relation dtu = — div(V^) and Jensen's inequality we get 

/ |[ div(VV') Illf-i dt = / HaHI^-i di > || / d t udt\\ 2 H -i = \\pi - po\\%-i- 
Jo Jo Jo 



Thus the Dacorogna-Moser transport is a minimum of the functional (|22p . 
so it is a transport minimizing the divergence of the velocity in //~ 1 -norm. 



4 Numerical Approximation of the 2D Optimal 
Extended Optical Flow 

The numerical method is based on a finite element time-space 1? least 
squares formulation (see (7j) of the transport problem (|6|). Define v n+1 
as 

and for a sufficiently regular function ip defined on Q, set 




dip dip dip V 
dt ' dx\ ' dx2 J ' 
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and 



div(^ +i 



dt^^dxi Kl 



Let • • • <£n} be a basis of a space-time finite element subspace 

V/i = {(f, piecewise regular polynomial functions, with 93(0, •) = (p(l, •) = 0}, 

for example, a brick Lagrange finite element of order one ([8]). Let Hh be the 
Lagrange interpolation operator. Let also Wh be the finite element subspace 
of Hq(Q), where the basis functions • ■ ■ iPm} are the traces at t = of 
basis functions {<fii}fL 1 . An approximation of problem is the following. 
For a discrete sequence of time t compute 



/ (^(V)(VM +1 -C"(t))|V4) dx 
Jn 



d tP Ut,-Whdx W h eW h , (24) 



and define v n+1 = Vp£ +1 . The L 2 least squares formulation of problem 
is defined in the following way. Consider the functional 



J(c h ) 



Q 



~n+l 



Vc h + v[n h ((i-t)p + t Pl )]) dxdt 



This functional is convex and coercive in the appropriate anisotropic Sobolev's 



space H = {cp e L 2 (Q); 



rn+l 



e L 2 (Q);v?(0,-) = cp(l, 



0} since 



rn+1 



L 2 (Q) 



IS 



the velocity field v n+1 is regular enough. Moreover 
a norm in H (see [7]). Set 

d=\i^\, e h = u h ({i-t)p + tp 1 ), 

and introduce 

K h = {<p h G Va; ^ + 0fc € I?} 

it is a closed convex subspace. Thus an approximation of ([6]) subject to the 
constraint p n+1 S D is defined with the following minimization problem: 

min J(c). (25) 

c£K h 



The minimizer of problem (|25p is 
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which is characterized by the Fermat's rule [2] 
v n+1 \Vc h ) ( v n+1 | Vip h ) dxdt> 



Q 



Q 



^_^+i|V^) (V +1 |V^) dx dt (26) 



for all iph € Tx h , the contingent cone to K^. Thus an approximation of the 
solution to problem © is 

pT x = c h + n fc ((i - t)po + t Pl ) e v h . 



Remark 4.1 Define the bilinear form a™ +1 (-, •) on Vh x oy 
a^+ 1 (( /3 ,^)= / (V +1 | VV> ) (V +1 | Vi> h ) dx dt. 



XTien £/te function Ch is the solution of 

as Zong as C/j is everywhere non negative. Moreover, if V h + consists of the 
non negative functions in Vh, the function Ch is the orthogonal projection 
of Qh onto V h + , according to the inner product induced by the bilinear form 

When the approximation with finite element of the algorithm proposed in 
section has converged, the computed solution ph can be decomposed as 
Ph = + P^+Oh, that is to say the approximated Dacoragna-Moser trans- 

port solution augmented with its orthogonal projection onto non negative 
functions space according to the inner product a v . 



In the next subsections, all the computations are done according to remark 

p n+ 

dn 



12.11 So it is assumed that the normal velocity satisfy v™ +1 = — = 0. 



4.1 The transport of a bump 

As a first example, the displacement of a bump is considered. More precisely, 
let U =] - l,l[x]-l,l[, y p = 0.5, r = 0.3, /3 > 0, and a > 0. For (x, y) G 
set r = sjx 2 + {y — yo) 2 , and define 

f p + ae-^/^o-^ if r 2 >r 2 
^° 1/3 elsewhere. 
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Then for r = \J x 2 + (y + yo) 2 , define 



Pi = 



{ 



P + ae" r2 /( r o- r2 ) if r 2 >r 2 



(28) 



(3 elsewhere, 



The domain O is subdivided into 40 x 40 elements, and the time interval 
]0, T[=]0, 1[ is subdivided into 60 elements, thus the linear system has 102541 
unknowns. Is is solved with a preconditioned conjugate gradient. 
In figures[T]l5l the shape of the bump is presented at the time-steps 0, 12, 24, 
36, 48, and 60, the first one is the initial shape, and the last one is the final 
shape. The next figure [T] describe the bump transport for (3 = 1. Remark 
that in this case this transport is very similar to the Dacorogna-Moser one. 
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Figure 1: Bump transport for = 1 



Then figures [21 [31 HI and [5] describe the bump transport for /3 = 0.5, 0.2, 
0.1, and 0.05. 
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Figure 2: Bump transport for = 0.5 
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Figure 3: Bump transport for = 0.2 
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Figure 4: Bump transport for f3 = 0.1 
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Figure 5: Bump transport for f3 = 0.05 



Remark that when j3 is small, the linear system becomes very ill-conditioned. 
Indeed for /3 < 0.1 usual preconditioners like the classical ICO one are useless. 
In this case the parallel Gram-Schmidt least squares preconditioner (DIAG 
+ LS CGS OPT) developed in [THJ [TTl [19] is used. Then the linear system is 
solved in parallel using 48 processors. The global tolerance for the iterative 
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scheme developed in section [2] is set to 0.01. 

The number of iterations, and the CPU time for each value of f3 is given in 
table I4TT1 For (3 < 0.05 the algorithm did not converge. 



p 


nb. iterations 


CPU time [s] 


residue 


1 


4 


128 


8.62 10" 4 


0.5 


5 


149 


9.23 10" 3 


0.2 


8 


201 


4.41 10" 3 


0.1 


10 


229 


9.47 10~ 3 


0.05 


55 


1019 


7.77 10~ 3 



Table 1: CPU time and number of iterations 



4.2 The left ventricle motion 

The iterative strategy described in Section 2 is then used to compute an 
approximated solution, and to reconstruct the systole to diastole images 
of a slice of a left ventricle. Ten time steps have been used to compute 
the solution, and 10000 degrees of freedom for the time-space least squares 
finite element. The approximated fixed point algorithm converges in about 
10 iterations with an accuracy of about 10~ 4 . In this case the usual ICO 
preconditioner is sufficient; this is essentially due to the fact that there is no 
large region in the domain with a very low density p. In the next figure [6J 
the initial image and the final image are presented. 




Figure 6: End of diastole of a left ventricular (a), of systole (b) 



In the following figure El two intermediate times 1/3 and 2/3 are shown. 
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Figure 7: Time step 3 and 6 



To summarize, in this work, we present a fixed point algorithm for the com- 
putation of the time dependent optimal mass transportation problem, al- 
lowing to handle the images tracking problem. The efficiency of the method 
has been tested with some 2D examples. 
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